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CFD  Studies  of  a  flapping  wine  NAV 


Introduction 


There  has  been  an  increased  level  of  activity  in  recent  years  in  the  development  of  small 
unmanned  vehicles.  Recent  interest  in  Micro  aerial  vehicles  (MAVs)  and  Nano  Air  Vehicles 
(NAVs)  with  wingspans  of  less  than  15  cm  and  7.5cm,  respectively,  has  focused  considerable 
attention  on  flapping  foil  propulsion.  In  the  low  Reynolds  number  (5000  <  Re  <  50,000)  regimes 
accessed  by  such  tiny  aircraft,  operating  at  very  low  speeds  (~10m/s),  flapping  wings  are 
generally  more  dynamically  advantageous  than  fixed  wing  configurations,  particularly  if 
operating  in  gusts  or  needing  to  avoid  obstacles.  Experimental  and  computational  studies  are 
also  providing  insight  into  the  aerodynamics  of  actual  biological  flight,  as  well  as  suggesting 
biomimetic  flapping  wing  configurations  that,  while  they  have  no  direct  counterpart  in  nature, 
may  be  suited  to  NAVs.  It  is  clear  that  a  large  number  of  creatures  achieve  their  desired  flight 
performance  characteristics  by  utilizing  either  active  or  passive  (or  both)  shape  modifications. 

This  research  effort  is  part  of  the  Phase  I  for  the  DARPA  program  on  NAVs.  The  goal  of 
this  effort  is  to  develop  an  NAV  with  a  wing  span  of  7.5cm.  The  NAV  will  be  subjected  to  gust 
environments  of  less  than  5  knots  and  majority  of  the  mission  will  be  in  hover.  In  the  Phase  I  of 
this  project,  multiple  models  and  prototypes  were  developed.  As  a  first  step  several  lightweight 
miniature  wings  and  linkages  were  created. 

The  computational  fluid  dynamics  study  focused  on  the  aerodynamic  efficiency  of  the  wings. 
This  included  a  significant  modeling  effort  in  order  to  study  the  complex  interaction  between  the 
wing  foils  and  the  wing  linkage.  The  majority  of  the  effort  was  directed  towards  gaining  insight 
to  the  design  of  high  performance  flapping  wings.  The  primary  focus  of  the  computational  study 
was  on  the  aerodynamic  performance  of  the  wings  in  order  to  obtain  a  near  optimal  kinematics  of 
the  wing  while  the  NAV  is  in  a  hover  mode.  The  secondary  effort  was  more  of  a  fundamental 
nature,  in  order  to  understand  the  differences  between  a  rotating  helicopter  blade  and  a  similarly 
designed  wing  undergoing  a  flapping  motion. 

An  unstructured  grid-based  unsteady  Navier-Stokes  solver  with  automatic  adaptive 
remeshing,  called  FEFLO,  was  used  to  compute  the  flow  about  the  flapping  wing  through  several 
complete  cycles  of  oscillation  for  each  of  the  cases  studied.  Ramamurti  and  Sandberg  [1] 
computed  2-D  flow  about  pitching  and  heaving  airfoils  and  compared  their  results  with  the 
experimental  work  of  Anderson  [2].  While  2-D  wing  section  investigations  can  yield  useful 
insights  on  the  coupled  pitching  and  heaving  dynamics,  it  is  essential  to  carry  out  3-D 
computations  to  understand  the  influence  of  spanwise  flow.  We  chose  to  pursue  the  fruit  fly  for 
computational  investigations  since  both  detailed  wing  kinematics  and  corresponding  force  data 
were  available.  Ramamurti  and  Sandberg  [3]  have  computed  the  force  production  in  Drosophila 
and  obtained  nearly  exact  agreement  with  the  measured  force  data  of  Dickinson  et  al.  [4] 
Ramamurti  and  Sandberg  [5]  have  also  carried  out  3-D  unsteady  computations  of  insect 
maneuvering.  The  computed  force  production  throughout  the  saccade  maneuver  again  was  in 
nearly  exact  agreement  with  measured  forces  of  Fry  et  al.  [6].  The  first  few  hover  cycles  of  this 
fruitfly  was  analyzed  and  yielded  a  figure  of  merit  of  20%. 


M  anuscript  approved  December  5,  2008. 
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Discussion  of  Results 


Version  1.0 


Initial  computations  were  performed  on  the  first  NAV  wing,  version  1.0,  designed  by  David 
Cylinder  of  Princeton  Plasma  Physics  Laboratory.  The  kinematics  from  the  RoboFly  experiment 
[4]  is  used  as  the  initial  kinematics.  The  wing  stroke  amplitude  was  scaled  by  a  factor  of  0.75 
from  the  original  RoboFly  kinematics,  Fig  lb.  The  wing  configuration  and  the  kinematics  are 
shown  in  Fig.  la-c.  3-D  unsteady  flow  simulations  (run  1.0)  were  conducted  at  a  flapping 
frequency  of  5Hz  and  80Hz  for  several  flapping  cycles.  At  a  flapping  frequency  of  5Hz,  the 
results  show  that  one  wing  produces  a  mean  lift  of  0.268  mN  and  a  mean  forward  force  of  nearly 
zero.  This  implies  that  two  wings  can  generate  enough  lift  to  hold  54mg.  As  the  frequency  is 
increased  to  80Hz,  Fig.  3,  the  mean  thrust  increases  by  a  factor  nearly  proportional  to  the  square 
of  the  frequency  ratios  to  0.0626  N  while  the  mean  forward  thrust  remaining  nearly  zero.  At  this 
frequency,  two  wings  can  produce  enough  lift  to  support  a  13g  vehicle.  The  time  history  of  the 
aerodynamic  power  input  to  the  wing  and  the  center  of  pressure  on  the  wing  were  computed  and 
are  shown  in  Fig.  4.  The  mean  power  input  to  one  wing  at  80Hz  is  nearly  1.01  W  and  the  mean 
center  of  pressure  is  located  at  70%  of  the  wing  span  from  the  hinge  point. 

Instead  of  the  version  1.0  wing  planform,  the  fruitfly  Drosophila  wing  shown  in  Fig.  5a  was 
used  (runl.l).  The  kinematics  were  maintained  to  be  the  same  as  in  Fig.l  and  results  showed 
that  there  was  a  slight  reduction  in  the  mean  thrust  produced  by  the  fruitfly  wing  to  0.0481  N, 
Fig  5c,  but  the  mean  aerodynamic  power  input  to  the  wing  is  reduced  by  nearly  40%  to  0.6  W, 
shown  in  Fig  5b.  Hence,  the  planform  was  modified  to  remove  a  small  region  near  the  trailing 
edge  portion  of  the  root  of  the  wing.  Figure  6  shows  a  few  modifications  to  this  wing. 


Fig.  1.  Wing  planform  and  kinematics  for  versionl.O. 
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Fig.  2.  Force  production  from  a  version  1.0  wing  flapping  at  5Hz  (runl.O),  (a)  time  history  of 
thrust  force,  (b)  thrust  for  2  flapping  cycles,  (c)  time  history  of  forward  force  and  (d) 
forward  force  for  2  flapping  cycles. 
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Fig.  3.  Force  production  from  a  version  1.0  (runl.O)  wing  flapping  at  80Hz. 
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Fig.  4.  (a)  Aerodynamic  power  input  and  (b)  center  of  pressure  on  a  wing  flapping  at  80Hz. 
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Fig.  5.  Results  of  a  flapping  Drosophila  wing,/=  80Hz,  (a)  planform  of  the  wing,  runl.l, 

(b)  aerodynamic  power  input  to  the  wing,  time  histories  of  (c)  thrust  and  (d)  forward  force. 
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NAV  2.0a  wing 


NAV  2.0c  wing  NAV  2.0d  wing 

Fig  6.  Modified  Wing  Shapes 

Kinematics  Extraction  from  video 


The  first  linkage  configuration  that  was  built  and  tested  was  a  double  crank  rocker  design. 
The  double  crank  rocker  was  a  planar  mechanism  and  had  the  advantage  of  being  simple,  which 
makes  it  amenable  for  miniaturization.  A  single  rotor  drove  a  four  bar  linkage  for  each  wing. 
However,  the  double  crank-rocker  mechanism  has  an  inherent  asymmetry  due  to  the  way  the  two 
four-bar  linkages  need  to  be  offset  when  attached  to  the  rotor. 

High  speed  video  was  obtained  for  this  mechanism  with  the  wing  planform  shown  in  Fig.  6a. 
One  of  the  slowest  runs  with  least  wing  twist  and  vibration  was  chosen.  Since  there  was  only 
one  camera  view  and  no  calibration  of  the  camera  was  available,  a  simple  digitization  of  the 
frames  was  performed.  On  each  frame  of  the  video  4  points  were  selected  as  shown  in  Fig.  7a. 
As  the  camera  view  was  normal  to  the  win  stroke  plane,  two  points  on  the  wing  root  (points  1 
and  2)  define  the  wing  stroke  angle.  The  wing  stroke  angle  (t)(t)  was  computed  directly  from  2-D 
pixel  coordinates  as  in  Eq.  la. 

In  order  to  obtain  the  angle  of  attack  or  the  wing  rotation  angle  two  points  were  selected  in 
the  mid  span  of  the  wing  (points  3  and  4),  one  at  the  leading  edge  and  the  other  on  the  trailing 
edge  of  the  wing.  The  angle  of  attack  was  computed  from  2-D  pixel  coordinates  using 
knowledge  of  the  pixel  distance  between  two  points  when  the  wing  plane  is  parallel  to  the  image 
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plane,  Fig.  7b,  and  by  scaling  the  distance  in  real  world  coordinates  (cm)  between  the  same 
points.  Using  the  video,  we  can  measure  the  pixel  distance  from  root  to  the  tip  of  the  wing.  We 
can  also  measure  the  distance  in  cm  from  root  to  the  tip  of  the  wing,  as  well  as  the  distance  in  cm 
from  point  3  to  point  4.  Knowing  the  distance  in  pixels  from  leading  to  trailing  edge  point  for 
each  frame  throughout  the  stroke,  we  can  calculate  the  angle  of  attack  a(t)  using  the  ratio  of 
measured  distance  and  the  distance  if  the  wing  is  parallel  to  the  view  plane,  as  in  Eq.  lb.  The 
wing  rotation  angle  measured  from  the  vertical  axis,  Y(t),  is  given  by  Eq.  Ic. 

The  kinematics  obtained  thus  is  shown  in  Fig.  8a  and  b.  The  3-D  unsteady  flow  was 
simulated  using  the  experimental  kinematics  shown  by  the  blue  curves  in  Fig.  8,  runl.2,  at  a 
wing  flapping  frequency  of  60Hz.  The  thrust  and  forward  force  time  history  is  shown  in  Fig.  9. 
These  traces  show  more  force  reversals  in  each  of  the  up  and  down  strokes.  Hence,  the  raw 
kinematics  obtained  from  the  video  is  smoothed  based  on  the  angular  accelerations.  The 
smoothed  kinematics  is  shown  by  the  green  curve  in  Fig.  8.  The  resultant  force  time  history  is 
shown  in  Fig.  10  and  we  can  see  the  high  frequency  oscillations  are  eliminated.  The  mean  thrust 
or  lift  per  wing  for  this  configuration  is  only  0.003  N,  shown  in  Fig.  1 1 .  Comparing  the  present 
kinematics  with  that  of  the  Robofly,  it  is  clear  that  the  wing  rotation  angle  is  not  properly  phased 
with  respect  to  stroke  reversal.  Hence,  the  kinematics  was  modified  so  that  the  Y(t)  closely 
matched  that  of  the  Robofly,  runl.2.2,  as  shown  by  the  black  curve  in  Fig.  8.  The  unsteady  flow 
simulation  using  this  modified  kinematics  showed  a  four-fold  improvement  in  the  lift  generated. 

Figure  1 1  shows  the  mean  lift  for  this  case  is  0.0131  N. 

/  \ 


^(l)  =  tan  '  — - - 


[fi-yi/ 


a(t)  =  sin  ' 


(la-c) 


Fig.  7.  Extracting  kinematics  from  the  video  of  the  double  crank-rocker  mechanism. 


Next,  a  laminar  viscous  solution  was  obtained  with  this  kinematics,  runl.2. 3,  to  study  the 
effect  of  viscosity.  The  Reynolds  number  for  this  case  is  approximately  9050,  based  on  the  mean 
wing  tip  velocity  and  a  mean  chord  of  1.73  cm.  It  is  clear  from  Fig.  11,  that  the  force  time  history 
remains  nearly  unchanged  for  this  case.  This  is  mainly  due  to  the  fact  that  the  forces  generated 
are  inertially  dominated  arising  from  the  large  motion  of  the  high  at  relatively  high  angles  of 
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attack.  Earlier  computations  on  the  Robofly  experiment  [3]  also  showed  similar  results  for  the 
effect  of  viscosity.  Hence,  in  this  study  inviscid  simulations  were  deemed  sufficient  for 
parametric  studies  and  optimization  of  the  wing  kinematics. 


Fig.  8.  Kinematics  obtained  from  the  double  crank  rocker  mechanism. 


Fig.  9.  Force  production  of  Version  1.0  with  experimental  kinematics, /=  60Hz,  runl.2. 
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time  (secs.) 


time  (secs.) 

Fig.  10.  Force  production  of  Version  1.0  with  smoothed  kinematics, /=  60Hz,  runl.2.1. 
Version  2.0 


The  next  generation  of  the  linkage  mechanism  that  was  developed  was  based  on  a  bicycle 
drive.  This  configuration  eliminated  the  asymmetries  present  in  the  double  crank  drive,  and  the 
right  and  left  wings  are  driven  symmetrically  with  in  manufacturing  tolerance,  while  keeping  the 
mechanism  is  robust  and  lightweight.  High  speed  video  of  the  wing  flapping  at  62. 3 Hz  were 
taken  by  Caltech  researchers  with  this  mechanism  for  wing  versions  b  and  c  shown  in  Fig.  6b 
and  c.  The  wing  and  the  linkage  mechanism  and  the  two  angles  of  rotations  are  shown  in  Fig. 
12a  and  b.  The  kinematics  was  obtained  by  digitizing  the  first  94  frames  of  NAV2.0c_v3.0  run 
and  is  shown  in  Fig.  13a  and  b.  The  kinematics  was  smoothed  based  on  the  angular  accelerations 
and  used  with  the  wing  planform  of  version  d  shown  in  Fig.  6d.  This  will  be  our  baseline 
configuration  to  improve  thrust  and  efficiency.  The  mean  thrust  for  one  wing  flapping  at 
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62.3Hz,  ran2.0  is  0.01 164  N,  shown  in  Fig.  14.  Note  that  there  is  only  one  peak  during  each  up 
and  downstroke. 


time  (secs.) 


time  (secs.) 

Fig.  11.  Force  production  of  a  wing  flapping  with  properly  phased  kinematics, /=  60  Hz. 

Parametric  Studies 


We  undertook  a  systematic  study  in  an  effort  to  the  study  the  effect  of  various  parameters 
that  govern  the  kinematics  in  an  attempt  to  optimize  the  kinematics  for  hovering  flight.  The 
possible  kinematic  parameters  include  stroke  angle  (t)(t),  the  wing  rotation  angle  Y(t),  the  wing 
flapping  frequency,  the  aspect  ratio  of  the  planform,  the  camber  and  twist  of  the  wing. 
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Fig.  12.  Bicycle  drive  mechanism,  NAV2.0c,  showing  the  two  wing  rotation  angles. 
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Fig.  13.  Kinematics  of  NAV2.0c. 

1 .  Effect  of  Stroke  angle  (c))): 


The  initial  stroke  angle  variation  (l)o(t)  is  obtained  from  NAV2.0c_v3.0  kinematics  (red  line  in 
Fig.  15).  The  amplitude  of  the  stroke  angle  is  varied  as  follows: 

=  1 .00  ^  =  104.6°  iNAV2.0c_v3.0) 

=  1. 25  =130.8°  ^  ^ 

=  1.45^0_-^^„  =151.7° 
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Fig.  14.  Thrast  time  history  for  the  baseline  configuration. 


Figure  16a  shows  the  variation  of  thrust  time  history  for  the  different  amplitudes  of  stroke  angle 
considered.  The  mean  thrust  increases  from  0.01 164N  to  0.0249N  by  increasing  the  stroke 
amplitude  from  the  baseline  value  of  104.6°  to  151.7°.  This  increased  amplitude  was  deemed 
mechanically  feasible  and  therefore  was  fixed  for  our  CFD  coupled  with  Design  of  Experiments 
simulations. 


The  thrust  produced  and  the  aerodynamic  power  were  computed  by  integrating  the  unsteady 
pressure  on  the  surface  of  the  wing  and  are  non-dimensionalized  as  follows: 

r  =  h  •  /  p  n  dA 


Paero=^-^  =  iP  W’ll  dA 


\pylK 

Cp  =  — ,  where 


(3a-e) 


where  p  is  the  computed  pressure  on  the  wing  surface,  w  is  the  velocity  of  the  wing,  n  is  the 
normal  vector  the  wing  surface,  h  is  the  unit  normal  in  the  hover  direction  and  Vtip  is  the  mean 
tip  speed  of  the  wing.  In  terms  of  the  non-dimensional  thrust  coefficient  all  the  three  curves 
nearly  collapse  into  one  as  shown  in  Fig.  16b.  The  mean  thrust  coefficient  per  flapping  cycle  is 
approximately  0.52. 
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Fig.  15.  Stroke  angle  variation. 
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Fig.  16.  Variation  of  thrust  and  eoefficient  of  thrust  time  history  with  stroke  amplitude. 
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1.  Effect  of  wing  rotation  angle  (y): 

The  baseline  wing  rotation  angle  variation,  Yo(t),  is  obtained  from  NAV2.0c_v3.0  kinematics 
(red  line  in  Fig.  17).  The  mid  stroke  wing  rotation  angle  was  varied  as  follows: 


y(0  =  y/„e*yo(0, 

=  1.00  =>  y_min  = -60°,  y_max  =  54°  {NAVlQc _v?)Q) 

=  0.75^y_min  =  -45°,  y_max  =  40.5°  ^  ’ 

=  0.50=>y_min  =  -30°,  y_max  =  27° 

Figure  18  show  the  variation  of  thrust  for  the  three  wing  rotations  angle  distributions  considered. 
As  the  mid  stroke  wing  rotation  angle  is  reduced  the  mean  thrust  increases  slightly  from  0.0249 
N  to  0.0253  N  mainly  due  to  the  improved  thrust  during  the  downstroke.  Further  reduction  in 
wing  rotation  angle  reduces  the  thrust  during  both  up  and  downstrokes  yielding  a  mean  value  of 
0.0207  N. 


0.000  0.005  0.010  0.015  0.020 

time  (secs.) 


Fig.  17.  Wing  rotation  angle  variation. 
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Fig.  18.  Thrust  variation  with  wing  rotation  angle,  run2.01  (red),  run2.04  (blue)  and  run2.07 
(green). 
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3.  Effect  of  Phasing  of  wing  rotation  angle  (|3): 

The  stroke  angle  for  this  study  will  be  with  (l)fac=  1.45.  The  NAV2.0e  kinematics  shown  in 
Fig.  19  (red  line  w.r.t.  black  line)  exhibits  delayed  rotation  in  the  wing  rotation  angle,  i.e.,  the 
wing  rotation  peaks  much  after  the  stroke  reversal.  Therefore,  the  wing  rotation  was  shifted  to  be 
in  phase  to  match  the  peak  of  the  stroke  angle  (blue  line).  This  was  termed  symmetric  with  a 
phase  angle  p  =  0°.  The  wing  rotation  was  further  modified  to  advance  the  rotation  (green  line), 
corresponding  to  a  p  =  15.3°.  The  base  line  Y(t)  with  y  fac  =  1.0  was  used  first. 

Figure  20  shows  the  thrust  variation  with  the  phasing  of  the  wing  rotation  angle  w.r.t.  the 
stroke  angle  and  the  effect  is  significant.  The  baseline  wing  rotation  (/off))  is  approximately 
15.3°  delayed  and  yields  a  mean  thrust  of  0.0249  N.  When  the  wing  rotation  is  symmetric,  the 
mean  thrust  improves  to  0.0284  N  and  advancing  the  wing  rotation  by  15.3°  improves  the  mean 
thrust  to  0.0294  N,  nearly  18%  improvement  over  the  baseline  wing  rotation.  In  the  advanced 
rotation  case,  the  appearance  of  a  second  thrust  peak  in  the  end  of  the  downstroke  or  beginning 
of  the  upstroke  suggests  the  wake  capture  effect. 
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Fig.  19.  Phasing  of  wing  rotation  angle  w.r.t.  stroke  angle. 
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Fig.  20.  Effect  of  phasing  of  wing  rotation  angle  on  the  thrust  time  history,  run  2.01  (red),  run2.02 
(blue)  and  run2.03  (green). 
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Further  parametric  studies  were  performed  varying  both  the  y  factor  and  the  phasing.  Table  1 
shows  the  variation  of  the  mean  thrust.  For  the  symmetric  and  advanced  cases,  reducing  the  mid 
stroke  wing  rotation  angle  results  in  a  reduction  in  the  mean  thrust.  These  results  show  that  we 
can  get  nearly  6gf  of  lift  from  2  wings  at  62.3Hz.  From  our  frequency  scaling  analysis,  we  can 
get  lOgf  at  80Hz.  These  results  were  used  by  the  Design  of  Experiments  (DOE)  analysis.  From 
that  analysis,  the  next  set  of  parametric  studies  increasing  the  wing  rotation  factor  to  1 .25  and  the 
phasing  to  30.6°  are  obtained. 


y  factor 

delayed  (-15.3°) 

Symmetric  (0°) 

advanced  (15.3°) 

1.00 

0.0249 

0.0284 

0.0294 

0.75 

0.0253 

0.0280 

0.0287 

0.50 

0.0207 

0.0224 

0.0219 

Table  1 .  Mean  Thrust  (N)  variation  with  mid  stroke  a  anc 

phasing  |3. 

4.  Effect  of  Frequency: 

In  order  to  test  the  effect  of  flapping  frequency  on  the  thrust/force  produced  and  the  power 
input  to  the  wing,  the  case  with  the  chord  factor  y  of  1.0  with  a  phase  shift  of  15.3°  and  a  rigid  fin 
was  selected.  The  frequency  was  increased  from  62.3Hz  to  72.5Hz  so  that  the  thrust  produced  by 
one  wing  is  approximately  5gf,  shown  in  Fig.  21. 


0.0000  0.0138  0.0276  0.0414  0.0552 

time  (secs.) 

Fig.  21.  Thrust  variation  at  a  flapping  frequency /=  72.5Hz.  run5.04.1. 

In  terms  of  the  non-dimensional  coefficients  the  thrust  and  the  power  distributions  for  the  two 
frequencies  tested  collapse  into  a  single  curve,  as  shown  in  Fig.  22a  and  b.  The  mean  thrust 
coefficient  is  approximately  0.61  and  the  mean  power  coefficient  is  0.83. 
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T 


time  (sec.) 

Fig.  22.  Variation  of  Thrust  and  Power  coefficients,  run5.04  (red)  and  run5.04.1  (blue). 

5.  Aspect  Ratio  (AR)  of  the  wing: 

For  this  study,  wing  version  d  shown  in  Fig.  6d  will  be  used.  The  AR  will  be  varied  by 
changing  the  max  chord.  The  baseline  max  chord  is  1.2cm.  This  will  be  changed  to  0.9cm  and 
1.5cm.  corresponding  to  a  chord  factor  of  0.75  and  1.25,  respectively.  The  phase  shift  for  this 
study  is  15.3°  and  the  y  factor  is  1.0.  The  results  shown  in  Fig.  23a  suggests  that  the  mean  lift 
varies  almost  linearly  with  the  chord  factor.  The  thrust  coefficients  for  these  three  cases  also 
merge  almost  into  a  single  curve  showing  the  linear  behavior  with  the  chord  factor. 
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Fig.  23.  Effect  of  chord  on  the  thrust  time  history,  run  2.03  (red),  run3.20  (green)  and  run3.21 
(blue). 

6.  Camber  and/or  Twist: 


For  reversing  camber,  an  approach  similar  to  our  earlier  study  of  the  BITE  vehicle  [7]  is 
used.  The  mean  camber  line  would  be  NACA  airfoil  camber.  The  maximum  camber  location  p 
is  fixed  to  be  0.25c  and  only  the  max  camber  ym.  Fig.  24a,  is  varied.  The  reversing  parameter 
A(t)  follows  the  rotation  angle  profile  of  the  Robofly  experiment,  shown  in  Fig.  lb.  The  wing 
shape  is  prescribed  throughout  the  cycle  at  several  control  points  at  the  leading  edge,  trailing 
edge  and  a  0.75  span  of  the  wing  from  the  root  location  as  shown  in  Fig.  24b. 
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y  .y^{2p.-F),  {0..  .p 

p  ^  ' 

- ^^(\-2p  +  2px  -x^\  {p<x  <\ 

c  Cj  c 

(5a-c) 

y{t)  =  A{t)-y^-C{z) 


A(0  =  - 


y  robotfly^^^ 


Similar  reversing  function  can  be  used  for  twist  with  0o  as  the  free  parameter  where  b  is  the 
maximum  span  of  the  wing  and  z  is  the  spanwise  location. 


d{z,t)  =  A{t)  20oz/h 


(6) 


For  the  case  with  the  y  factor  of  1 . 125,  and  phase  shift  of  7.65°,  a  maximum  camber  of  10% 
chord  produced  a  higher  mean  thrust  of  0.038  N  compared  to  the  rigid  wing  which  produced 
approximately  0.0339N,  as  shown  in  Fig.  25. 


(b) 

Mp  =  0.25,  ym  =  0.05 

■  p=  0.25,  ym  =  0.01 

■  rigid  body  rotation 


Fig.  24.  Prescribing  reversing  camber  kinematics,  (a)  wing  cross  section  and  (b)  control  points. 
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time  (secs.) 

Fig.  25.  Effect  of  reversing  camber  on  the  thrust  production. 

7.  Optimization  using  CFD  and  Design  of  Experiments: 

For  the  optimization  of  the  various  parameters  that  control  the  kinematics  of  the  wing,  the 
kinematics  obtained  from  the  experimental  observation  was  used  as  a  baseline.  We  then  varied 
the  other  parameters  using  scaling  factors  from  the  baseline  kinematics,  such  as  the  wing  stroke 
amplitude,  4>fac,  the  wing  rotation  angle  or  mid-stroke  angle  of  attack,  y/ac,  the  phasing  between 
the  wing  rotation  and  stroke  reversal,  13,  and  the  maximum  chord  length  and  the  maximum 
camber  at  %  chord  location,  ym- 

During  the  initial  stages  of  the  optimization  study,  it  was  found  that  the  thrust  generated 
increases  as  the  stroke  amplitude  was  increased  from  the  baseline  value  to  151.7°  which  was 
deemed  the  maximum  feasible  for  the  current  mechanism.  This  value  corresponds  to  a  (pfac  of 
1.45  and  was  fixed  for  the  remainder  of  the  simulations.  Next,  from  our  frequency  scaling 
computational  experiments,  it  was  found  that  the  thrust  produced  T,  scales  proportional  to and 
the  aerodynamic  power  input,  Paem,  scales  proportional  to /^.  Hence,  the  flapping  frequency  was 
fixed  at  the  baseline  value  of  62.3  Hz.  The  maximum  thrust  of  0.0396  N  was  obtained,  run  5.10, 
with  a  Yfac  =  1.05925,  (3  =  22.95°  and  =  0.1.  The  figure  of  merit  at  this  condition  was  only 
18%. 


The  objective  of  the  optimization  is  to  maximize  the  figure  of  merit  (FOM)  while  producing 
a  10  gf  of  lift.  In  order  to  compare  FOMs  between  simulations,  first  the  frequency, /yo,  that 
would  produce  10  gf  of  lift  was  found  using  the  computed  thrust  T  and  the  scaling  laws  for  the 
frequency.  Next,  the  aerodynamic  power  Paero  is  scaled  to  Pareo, constrained  which  is  the  power 
required  to  produce  10  gf  of  lift.  Then,  the  FOM  is  obtained  from  this  power  and  the  ideal 
induced  power  from  an  equivalent  actuator  disk.  This  is  shown  in  the  following  equation: 
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/lo  =  f^mg/2T 


P  =2P 

aero  constrained  aero 


=  p  IP 

induced  aero  constrained 


(7) 


(mgf^ 

^induced  ^  A  \l/2 

swept) 

2  3 

where  m  is  the  mass  of  the  vehicle  (=  0.01  kg),  g  =  9.8  m/s  ,  p  is  the  density  of  air  (=1.2  kg/m  ) 
and  Aswept  is  the  area  swept  by  the  wing  depending  on  the  stroke  amplitude. 


The  set  of  parameters  for  each  computational  simulation  was  chosen  using  statistical  Design 
of  Computer  Experiments  [8]  method  to  find  the  optimal  configuration  while  minimizing  the 
number  of  simulations.  From  the  computational  experiments,  it  was  found  that  the  near  optimal 
set  of  kinematics  parameters  is  jfac  =  1.25,  yS  =  -15.3°  and  gm  =  0.04.  From  the  computational 
results  varying  the  set  of  the  three  control  parameters  and  using  the  design  of  experiments,  it  was 
found  that  the  maximum  FOM  can  be  obtained  at  y/ac  =  1.256, 13  =  -1 1.02°  andgm  =  0.0455,  and 
the  maximum  FOM  was  36%.  Figure  26  shows  the  contours  of  FOM  at  3  different  camber 
depths.  The  optimal  design  point  is  shown  by  the  star  in  Fig.  24b,  for  convenience,  although  it 
occurs  on  a  different  plane  in  the  3-parameter  search  space. 


For  the  near  optimal  configuration,  the  FOMkover  was  found  to  be  34.5%.  For  this  case,  the 
mean  thrust  generated  by  one  wing  is  0.0241  N  and  the  Paero  is  0.1649  W.  The  thrust  and  power 
time  histories  are  shown  in  Fig.  27  a  and  b,  respectively.  Figure  27a  also  shows  that  a  rigid  wing 
{ym  =  0.0),  produces  a  reduced  mean  thrust  of  0.02196  N  yielding  an  FOM  of  29.4%.  Figure  28 
shows  the  pressure  distribution  on  the  top  and  bottom  surfaces  of  the  wing  corresponding  to  the 
instant,  t  =  0.0  at  which  maximum  thrust  is  obtained.  It  is  clear  that  lower  suction  peaks  and 
larger  high  pressure  peaks  on  the  bottom  and  top  surfaces,  respectively,  are  attained  in  the 


0.625  Yfac  1.875 


(a)  ym  =  0.0  (b)  y^  =  0.05  (c)  =  0. 1 

Fig.  26.  Contours  of  FOM  using  Design  of  Experiments  to  find  optimal  combinations  of  design 
parameters  of  the  kinematics  of  the  flapping  wing. 
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Fig.  27.  Variation  of  (a)  thrust  and  (b)  aerodynamic  power  for  the  near  optimal  configuration, 
run6.05  (red)  and  run6.06  (blue). 
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Fig.  28.  Surface  pressure  distribution  on  the  rigid  and  flexible  wings  at  t  = 
upstroke. 
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1  St  Cycle  2nd  Cycle 

Fig.  29.  Surface  pressure  distribution  during  the  first  and  second  flapping  cycles  at  the  instant  of 
maximum  thrust  during  the  down  stroke,  showing  the  effect  of  wake  capture. 
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Fig.  30.  Pressure  and  velocity  vectors  on  a  plane  cut  z  =  3.5  cm. 
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flexible  wing  case  compared  to  the  rigid  case,  resulting  in  a  larger  lift  force.  Also,  it  is 
interesting  to  note  that  the  peak  thrust  during  the  first  cycle  of  flapping,  t  =  0.003s  is  larger 
compared  to  the  peaks  in  the  successive  cycles.  This  is  mainly  due  to  the  influence  of  the 
downwash  in  the  wake  which  reduces  the  effective  angle  of  attack.  Fig.  29.  This  effect  is  also 
observed  in  experiments  by  Birch  and  Dickinson  [9].  Figure  30  a-d  show  the  pressure  and 
velocity  vectors  on  a  plane  3.5cm  from  the  symmetry  plane  at  these  instants  during  the  first  two 
flapping  cycles.  It  is  clear  that  the  leading  edge  vortex  during  the  first  cycle  is  much  stronger 
compared  to  the  second  cycle.  Figure  31  show  the  particles  traces  and  the  presence  of  the 
leading  edge  vortex  which  merge  with  the  wing  tip  vortex  close  to  the  instant  when  the 
maximum  thrust  is  attained  during  the  upstroke. 


Fig.  31.  Particle  traces  at  t  =  0.0298  s. 

8.  Effect  of  Planform: 


Several  plan  forms  of  the  wing  have  been  tested  both  experimentally  and  in  computational 
studies,  as  shown  in  Fig.  32.  The  characteristic  dimension  that  distinguishes  these  planforms  is 
the  mean  chord  of  the  wing  which  can  be  defined  as 

C  =  (8) 

L 

span 

where  is  the  area  of  the  wing  and  Lspan  is  the  span  of  the  wing.  In  all  our  studies,  the  span  of 
the  wing  is  maintained  to  be  3.75cm. 

The  initial  configuration  of  the  wing  Fig.  32b  used  with  the  first  mechanism  built  is  found  to 
be  restrictive  in  the  stroke  amplitude  and  therefore,  the  trialing  edge  portion  near  the  root  was 
modified.  Version  2.0c,  Fig.  32c,  is  further  modified  to  remove  the  leading  edge  portion  of  the 
wing  ahead  of  the  wing  rotation  axis  to  yield  better  turning  characteristics  during  stroke  reversal, 
resulting  in  wing  version  d.  Fig.  32d.  The  optimization  study  for  the  kinematics  is  performed 
with  the  version  d  of  the  wing  and  the  near  optimal  kinematics  is  used  to  evaluate  the 
performance  of  the  two  paddle-shaped  wings  Fig.  32e  and  f  Finally,  the  planform  of  the  fruitfly 
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wing  is  scaled  to  have  the  3.75cm  wing  span  and  is  used  to  compare  the  performance  of  the 
version  d  of  the  wing  using  the  kinematics  obtained  from  the  Robofly  experiments. 


Fig.  32.  Various  planforms  of  wings  tested  in  experiments  and  computations. 

(a)  fruitfly  Drosophila,  C  =  0.8524  cm,  (b)  versionl.O,  C  =  1.4291  cm,  (c)  version  2.0, 
C  =  1.509  cm,  (d)  version  d,  C  =  0.6944  cm,  (e)  paddle  1,  C  =  0.7619  cm  and 
(f)  paddle  2,  C  =  0.8755  cm. 

With  the  near  optimal  kinematics  described  in  the  previous  section,  the  wing  planform  was 
changed  to  the  paddle  2,  Fig.  32f.  This  resulted  in  a  mean  trust  of  0.03048  N,  Fig.  33a,  and  a 
mean  aerodynamic  power  of  0.23  W,  Fig.  33b,  yielding  a  slightly  improved  FOM  of  35.15%. 

Propeller  vs  Flapping 


In  order  to  study  the  differences  between  the  propulsive  efficiency  of  a  propeller  to  a 
flapping  wing,  first  a  propeller  was  designed  by  Steve  Morris  to  produce  lOg  of  thrust.  This 
propeller  has  a  high  twist  from  the  root  to  the  tip.  The  shape  of  the  propeller  is  given  as 


c  =  0.75  -z/3  (in) 

^  ,  _i/0.23873\ 

0„.„,(z)  =  tan  M — - — I 


prop^ 


-3.0  (deg),  0< 


where  c  is  the  chord,  z  is  the  radius  of  the  spanwise  location  and 
plane  of  rotation.  The  maximum  camber  at  each  spanwise  section 
25%  of  the  chord.  The  resulting  geometry  is  shown  in  Fig.  34. 


z^l.5". 


(9) 


dprop  is  the  twist  angle  from  the 
is  set  to  4%  of  the  chord  and  at 
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Fig.  33.  Thrust  and  aerodynamic  power  variation  at  the  near  optimal  kinematics  using  paddle  2 
planform,  run  7.09 


Leading  Edgej 


j  Trailing  edge 

root 


(0  =  87°) 

Fig.  34.  Three  views  of  the  geometry  of  the  rotating  propeller. 
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Computations  were  performed  for  this  propeller  rotating  at  8750  rpm  in  a  rotating  frame  of 
reference.  The  resultant  thrust  generated  by  the  two  blades  is  0.10817  N  and  the  aerodynamic 
power  input  to  the  propeller  is  0.6405  W,  resulting  in  a  FOM  of  approximately  54%.  The 
pressure  and  the  magnitude  of  velocity  on  the  top  and  bottom  surfaces  of  the  propeller  blade  are 
shown  in  Fig.  35.  The  same  configuration  was  computed  in  a  stationary  frame  of  reference  with 
the  propeller  blades  modeled  as  moving  bodies  at  the  specified  rpm.  Figure  36  shows  the  time 
history  of  the  thrust  and  aerodynamic  power  input  for  this  propeller.  The  mean  thrust  for  the  2 
bladed  propeller  is  0.16761  N  and  the  mean  power  input  is  1.0053  W,  resulting  in  a  FOM  of 
66.3%.  Again,  this  is  similar  to  the  results  obtained  with  the  helicopter  blade  design  code  of 
Steve  Morris. 


a.  pressure  contours 


b.  magnitude  of  velocity 

Fig.  35.  Surface  pressure  and  velocity  distribution  on  the  propeller  blades  rotating  at  8750  rpm. 

In  order  to  study  the  effect  of  viscosity,  a  viscous  computation  was  carried  out  at  a  Re  = 
34,000  based  on  the  mean  tip  speed  of  34.36m/s  and  the  mean  chord  of  1.479cm.  Figure  37 
shows  the  time  history  of  thrust  and  the  aerodynamic  power  input  to  the  propeller.  The  overall 
figure  of  merit  for  this  case  is  reduced  to  59.6%  from  the  inviscid  value  of  66.3%. 

Propeller  as  a  flapper 

The  propeller  blade  described  above  was  used  as  a  flapping  wing  with  reversing  camber  and 
twist.  The  robot  fly  rotation  angle  profile  (Eq.  6c)  is  used  to  provide  the  reversing  camber  and 
twist. 

d{z,t)  =  A{t)  (10) 
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Fig.  36.  Thrust  and  aerodynamic  power  variation  in  a  propeller  rotating  at  8750  rpm. 


The  kinematics  for  the  flapping  propeller  blade  is  shown  in  Fig.  38.  For  this  blade  flapping  at 
145.8333  Hz,  the  resultant  time  history,  Fig.  39a,  shows  that  the  mean  thrust  produced  by  one 
wing  is  0.171  N  and  the  mean  aerodynamic  power  input  to  the  blade  is  6.2776  W,  Fig.  39b, 
resulting  in  a  FOM  of  approximately  15.5%.  Although  the  mean  thrust  produced  by  the  flapping 
propeller  blade  is  almost  twice  that  of  the  rotating  propeller  (=  0.0838  N),  the  mean  power 
required  per  cycle  is  nearly  12  times  that  of  the  propeller  (=  0.5027  W)  per  blade,  resulting  in  a 
reduced  FOM.  From  the  frequency  scaling,  we  can  obtain  a  thrust  of  10  gf  with  two  wings 
flapping  at  78  Hz.  Figure  39a  shows  that  the  thrust  production  is  negative  just  after  stroke 
reversals  and  from  Fig.  39b  we  can  see  that  the  power  input  reaches  a  maximum  just  before 
stroke  reversals.  Similar  trends  are  shown  in  the  time  history  of  the  coefficient  of  thrust  and  drag 
shown  in  Fig.  40.  The  ratio  of  the  thrust  and  drag  coefficients.  Fig.  41,  show  that  after  stroke 
reversals  the  ratio  reaches  a  large  negative  values.  At  these  instants  the  thrust  produced  is 
negative.  This  is  largely  due  to  the  specified  kinematics,  which  shows  a  large  acceleration  of  the 
trailing  edge  points  in  the  y-direction,  Fig.  42.  This  produces  a  high  pressure  on  the  top  surface 
of  the  blade  and  a  lower  pressure  on  the  bottom,  resulting  in  the  negative  thrust/lift. 
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Fig.  37.  Thrast  and  aerodynamic  power  variation  in  a  propeller  rotating  at  8750  rpm,  Re  ^ 
34,000. 


Fig.  38.  Kinematics  of  the  flapping  propeller,  ran7.10. 
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Fig.  39.  Thrust  production  and  aerodynamic  power  input  to  the  flapping  blade, /=  145.8333Hz. 


The  flow  around  the  flapping  foil  at  instants  just  before  and  after  stroke  reversals  and  in  the 
middle  of  the  upstroke  was  examined.  The  pressure  and  velocity  distributions  on  two  spanwise 
plane  cuts  were  obtained.  Figure  43  shows  the  locations  of  the  plane  euts  at  1/3  and  2/3  along  the 
span  from  the  root  and  the  rake  with  respeet  to  the  blade.  Fig.  44a  and  b  show  the  pressure 
distribution  and  Figs.  44c  and  d  show  the  magnitude  of  veloeity  on  these  plane  cuts  at  the  instant 
when  the  thrust  generation  is  maximum  right  before  stroke  reversal  at  t  =  0.0046  s.  It  is  elear 
from  Fig.  44d  and  e,  that  both  leading  and  trailing  edge  vortiees  are  present.  Figure  45  shows  the 
partiele  traces  released  from  the  rake  and  the  presence  of  the  leading  edge  vortex  and  the  wing 
tip  vortex  are  seen.  The  pressure  distribution  on  the  top  and  bottom  surfaces  of  the  wing  shows 
high  pressure  on  the  bottom,  Fig.  45c,  and  lower  pressure  on  the  top.  Fig.  45d,  leading  to  a  peak 
thrust  foree. 
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Fig.  40.  Time  history  of  the  coefficient  of  thrust  and  drag. 
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Fig.  41.  Time  history  of  Ct/Cd  for  the  flapping  propeller  blade. 
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Fig.  42.  Vertical  acceleration  of  2  trailing  edge  points. 

With  the  flapping  frequency  reduced  to  80Hz,  the  mean  thrust  developed  by  this  propeller 
blade  as  a  flapper  is  0.04725N  and  the  mean  aerodynamic  power  is  0.971  W,  as  shown  in  Fig. 
46a  and  b,  respectively.  It  is  clear  that  the  aerodynamic  power  peaks  near  stroke  reversals. 
Therefore,  the  mean  aerodynamic  power  and  the  mean  thrust  generated  during  the  middle  30% 
of  each  of  the  up  and  down  strokes  were  computed  and  the  resulting  FOM  is  48%.  Hence,  the 
kinematics  was  modified  so  that  the  stroke  reversal  for  twist  is  achieved  faster  as  shown  in  Fig. 
47.  Also  for  this  case,  the  stroke  amplitude  was  reduced  to  160deg  and  the  chord  at  the  root  of 
the  blade  was  reduced  by  20%.  These  modifications  were  necessary  in  order  to  avoid  the  trailing 
edges  of  the  blade  from  colliding  into  each  other  during  stroke  reversals.  With  this  set  of 
kinematics,  the  mean  thrust  from  the  blade  flapping  at  80.0Hz  is  0.0441N  and  the  mean  power 
input  is  0.563W,  resulting  in  an  overall  FOM  of  nearly  24%;  the  FOM  for  the  middle  80%  of  the 
stroke  is  66.7%  which  matches  that  obtained  from  the  rotating  propeller.  The  thrust  and  power 
time  histories  for  this  case  are  shown  in  Fig.  48a  and  b.  Attempts  to  change  the  plan  form  of  the 
blade  to  a  rectangular  one  and  modifying  the  phasing  of  the  twist  with  respect  to  stroke  reversal 
by  advancing  and  delaying  by  8%  of  the  stroke  period  did  not  result  in  any  improvement  in  the 
FOM.  Finally,  the  stroke  profile  was  modified  from  a  triangular  waveform  profile  to  a  smoothed 
profile  that  is  obtained  from  the  Robofly,  as  shown  in  Fig.  49a.  The  twist  and  the  camber  for  this 
case  are  as  shown  in  Fig.  50b  and  c.  With  this  modification,  the  mean  thrust  over  a  cycle  is 
0.03427N  and  the  mean  power  input  is  0.03437  W,  yielding  an  overall  FOM  of  nearly  27%.  The 
corresponding  FOM  for  the  middle  60%  of  the  stroke  is  57%.  The  thrust  and  power  time 
histories  for  this  case  are  shown  in  Figs.  50a  and  b,  respectively.  Another  attempt  is  made  to 
improve  by  changing  the  mid  stroke  angle  of  attack  at  the  tip  of  the  blade,  see  Fig.  34,  from  6°  to 
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15°.  The  mean  thrast  in  this  case  is  0.0433 1  N  and  the  mean  input  power  is  0.04758  W,  resulting 
in  an  overall  FOM  of  27.6%,  but  the  mid-stroke  FOM  dropped  to  44.5%. 


(b) 


Fig.  43.  Plane  cuts  and  particle  rake  locations  for  flow  visualization. 


The  effect  of  viscosity  was  studied  for  the  flapping  propeller  case,  in  order  to  examine  if 
there  is  separated  flow  along  the  chord  during  the  mid  stroke.  The  kinematics  for  this  case  is 
same  as  run7.20,  shown  in  Fig.  49.  The  Reynolds  number  based  on  the  wing  tip  velocity  and  the 
men  chord  of  1.09  cm  is  12,200.  The  mean  thrust  and  power  input  are  shown  in  Fig.  51.  The 
mean  thrust  dropped  slightly  from  0.03029  N  and  the  mean  power  per  cycle  also  reduced  to 
0.3139  W,  resulting  in  an  FOM  of  24.8%;  the  mid  stroke  FOM  for  this  case  is  also  slightly 
reduced  to  53.9%. 
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Fig.  44.  Pressure,  magnitude  of  velocity  and  velocity  vectors  at  t  =  0.0046  secs. 


A  closer  look  at  the  flow  shows  that  near  the  stroke  reversals  both  leading  and  trailing  edge 
vortices  are  being  shed  at  stroke  reversals,  t  =  0.0134s  and  t=0.0196s,  and  during  midstroke,  t  = 
0.0165s  and  t  =  0.0216s,  only  tip  and  root  vortices  are  present  around  the  blade.  These  features 
are  shown  in  the  iso-vorticity  contours  in  Fig.  52  a-d.  Figure  53  shows  the  second  invariant  of 
the  gradient  of  velocity,  a  quantity  that  is  a  good  indicator  of  the  vortex  core  [10].  From  this 
figure  also,  we  see  the  leading  and  trailing  edge  vortices  and  the  tip  vortex. 
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Fig.  45.  Flow  past  a  flapping  propeller  blade  at  t  =  0.0046  secs,  (a)  position  of  the  blade  during 
the  downstroke,  (b)  particle  traces  near  the  leading  edge,  surface  pressure  distribution  on 
(c)  top  and  (d)  bottom  of  the  blade. 

Summary  and  Future  Work 


Computational  fluid  dynamics  studies  of  the  3-D  unsteady  flow  past  a  flapping  wing  were 
conducted  for  the  development  of  a  NAV.  The  focus  of  this  computational  study  was  to  gain 
insight  to  the  design  of  high  performance  flapping  wings.  The  primary  focus  of  the 
computational  study  was  on  the  aerodynamic  performance  of  the  wings  in  order  to  obtain  a  near 
optimal  kinematics  of  the  wing  while  the  NAV  is  in  a  hover  mode.  CFD  simulations  coupled 
with  Design  of  Experiments  were  conducted  to  arrive  at  the  near  optimal  kinematics  for  the 
NAV  design.  Several  parametric  studies  were  performed  varying  the  stroke  angle,  the  wing 
rotation  angle  or  the  angle  of  attack,  the  phasing  between  these  two  rotations,  the  flapping 
frequency,  the  aspect  ratio  of  the  wing  and  the  planform.  It  was  found  that  the  near  optimal 
kinematics  yielded  an  efficiency  or  a  figure  of  merit  in  hover  of  35%.  In  comparison,  from  our 
previous  work  on  hovering  fruitflies,  the  FOM  in  hover  is  only  around  20%.  The  lower 
efficiency  in  flapping  propulsion  may  be  due  to  the  losses  generated  at  stroke  reversals. 
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Fig.  46.  Thrast  production  and  aerodynamic  power  input  to  the  flapping  blade, /=  80.0Hz. 


Fig.  47.  Improved  reversing  twist  kinematics  for  a  flapping  propeller  blade,  ran7.14. 
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Fig.  48.  Thrast  production  and  aerodynamic  power  input  to  the  flapping  blade, /=  80.0Hz. 


Fig.  49.  Modified  stroke  angle  profile  for  the  flapping  propeller  blade,  ran7.20. 
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Fig.  50.  Thrast  production  and  aerodynamic  power  input  to  the  flapping  blade  with  modified 
stroke, /=  80.0Hz. 

The  secondary  effort  was  more  of  a  fundamental  nature,  in  order  to  understand  the 
differences  between  a  rotating  helicopter  blade  and  a  similarly  designed  wing  undergoing  a 
flapping  motion.  We  found  that  an  optimally  designed  propeller  blade  was  able  to  operate  with 
an  FOM  of  60%  while  the  similarly  shaped  wing  flapping  at  the  same  frequeney  attained  only  an 
overall  FOM  of  25%.  Nevertheless,  the  flapping  wing  at  the  same  frequency  as  the  rotating 
propeller  was  able  to  generate  70%  higher  thrast.  Again,  the  lower  efficieney  values  are  due  to 
the  eontinuous  acceleration  and  deeeleration  of  the  flow.  Also,  the  flapping  propeller  blade 
attains  a  similar  FOM  of  60%  during  the  midstroke  of  the  cycle.  This  suggests  that  the  propeller 
blade  which  was  optimized  for  the  midstroke  angle  of  attaek  was  not  well  suited  at  the  stroke 
reversals,  during  whieh  the  wing  undergoes  large  ehanges  in  the  angle  of  attack.  The  question  of 
whether  or  not  a  flapping  wing  in  hover  is  neeessarily  less  efficient  than  a  helicopter  propeller  is 
still  open  due  to  the  fact  that  certain  kinematic  patterns  have  to  be  pieked.  It  is  not  elear  if  poorer 
performance  is  fundamental  to  flapping  or  an  artifact  of  the  ehosen  kinematics.  Further  studies 
are  needed  to  eliminate  the  adverse  effects  of  vorticity  shedding  seen  at  stroke  reversals. 
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Another  study  that  is  pertinent  to  hovering  vehicles  with  flapping  wings,  is  to  employ  a  four¬ 
winged  configuration  similar  to  that  of  a  Dragonfly.  This  opens  up  an  area  to  investigate  the 
effect  of  the  phasing  of  the  hind  wing  with  respect  to  the  forewing  in  order  to  augment  the 
performance. 
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Fig.  5 1 .  Thrust  production  and  aerodynamic  power  input  to  the  flapping  blade,  Re  =  12,200. 
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Fig.  52.  Iso-vortieity  contours  during  a  flapping  cycle, /=  80  Hz. 
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(a)  t  =  0.0134s 


WTV 


(b)t  =  0.0165s 


(c)t  =  0.0179s  (d)t  =  0.0196  s 

Fig.  53.  Vortical  activity  in  the  flow  depicted  by  the  leading  edge  vortex  (LEV),  trailing  edge 

vortex  (TEV)  and  wing  tip  vortiees  (WTV).  Iso-eontours  of  the  seeond  invariant  of  Vv 
are  shown. 
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Table  2.  Physical  and  Kinematic  parameters  for  CFD  simulations 


run  # 

Wing 
radius  (m) 

Radius  of 
non-wash 
region  (m) 

Flapping 
Freq  (Hz) 

(l)max  - 
(l)min  (deg) 

Chord 

Factor 

yfactor 

Phase 
Shift:  |3 
(deg) 

Max. 

camber 

1.01 

0.0425 

0.0115 

80.0 

120.0 

1.66667 

0.00 

0.000 

1.1 

0.0425 

0.0115 

80.0 

120.0 

1.30000 

0.00 

0.000 

1.2.2 

0.0369 

0.0042 

60.0 

101.0 

1.66667 

0.00 

0.000 

1.2.3 

0.0369 

0.0042 

60.0 

101.0 

1.66667 

0.00 

0.000 

2.01 

0.0375 

0.0065 

62.3 

151.7 

1.00000 

1.00000 

-15.30 

0.000 

2.02 

0.0375 

0.0065 

62.3 

151.7 

1.00000 

1.00000 

0.00 

0.000 

2.03 

0.0375 

0.0065 

62.3 

151.7 

1.00000 

1.00000 

15.30 

0.000 

2.04 

0.0375 

0.0065 

62.3 

151.7 

1.00000 

0.75000 

-15.30 

0.000 

2.05 

0.0375 

0.0065 

62.3 

151.7 

1.00000 

0.75000 

0.00 

0.000 

2.06 

0.0375 

0.0065 

62.3 

151.7 

1.00000 

0.75000 

15.30 

0.000 

2.07 

0.0375 

0.0065 

62.3 

151.7 

1.00000 

0.50000 

-15.30 

0.000 

2.08 

0.0375 

0.0065 

62.3 

151.7 

1.00000 

0.50000 

0.00 

0.000 

2.09 

0.0375 

0.0065 

62.3 

151.7 

1.00000 

0.50000 

15.30 

0.000 

5.01 

0.0375 

0.0065 

62.3 

151.7 

1.25000 

0.93745 

15.30 

4.485 

5.02 

0.0375 

0.0065 

62.3 

151.7 

1.25000 

1.12500 

-7.65 

0.000 

5.03 

0.0375 

0.0065 

62.3 

151.7 

1.25000 

0.62500 

22.95 

0.000 

5.04 

0.0375 

0.0065 

62.3 

151.7 

1.25000 

1.00000 

15.30 

0.000 

5.05 

0.0375 

0.0065 

62.3 

151.7 

1.25000 

1.25000 

22.95 

0.000 

5.06 

0.0375 

0.0065 

62.3 

151.7 

1.25000 

0.62500 

-7.65 

0.402 

5.07 

0.0375 

0.0065 

62.3 

151.7 

1.25000 

0.62500 

22.95 

10.000 

5.08 

0.0375 

0.0065 

62.3 

151.7 

1.25000 

0.89846 

22.95 

4.993 

5.09 

0.0375 

0.0065 

62.3 

151.7 

1.25000 

1.12500 

7.65 

10.000 

5.10 

0.0375 

0.0065 

62.3 

151.7 

1.25000 

1.05925 

22.95 

10.000 

5.11 

0.0375 

0.0065 

62.3 

151.7 

1.25000 

1.12500 

-7.65 

5.829 

5.12 

0.0375 

0.0065 

62.3 

151.7 

1.25000 

0.62500 

7.65 

5.390 

5.13 

0.0375 

0.0065 

62.3 

151.7 

1.25000 

0.89797 

7.65 

0.000 

5.14 

0.0375 

0.0065 

62.3 

151.7 

1.25000 

0.74633 

-7.65 

10.000 

5.15 

0.0375 

0.0065 

62.3 

151.7 

1.25000 

0.86646 

7.65 

5.019 

5.09.1 

0.0375 

0.0065 

62.3 

151.7 

1.25000 

1.12500 

7.65 

0.000 

5.10.1 

0.0375 

0.0065 

62.3 

151.7 

1.25000 

1.05925 

22.95 

0.000 

5.04.1 

0.0375 

0.0065 

72.5 

151.7 

1.25000 

1.00000 

15.30 

0.000 

5.04.2 

0.0375 

0.0065 

62.3 

151.7 

1.22250 

1.00000 

15.30 

0.000 

5.04.3 

0.0375 

0.0065 

62.3 

151.7 

1.66667 

1.00000 

15.30 

0.000 

5.04.4 

0.0450 

0.0065 

62.3 

151.7 

1.55000 

1.00000 

15.30 

0.000 

6.01 

0.0375 

0.0065 

62.3 

151.7 

1.25000 

1.50000 

15.30 

0.000 

6.02 

0.0375 

0.0065 

62.3 

151.7 

1.25000 

1.50000 

0.00 

0.000 

6.03 

0.0375 

0.0065 

62.3 

151.7 

1.25000 

1.50000 

-15.30 

0.000 

6.04 

0.0375 

0.0065 

62.3 

151.7 

1.25000 

1.50000 

-30.60 

0.000 

6.05 

0.0375 

0.0065 

62.3 

151.7 

1.25000 

1.25000 

-15.30 

0.000 

6.06 

0.0375 

0.0065 

62.3 

151.7 

1.25000 

1.25000 

-15.30 

4.000 

7.01 

0.0375 

0.0065 

62.3 

120.0 

1.25000 

0.00 

0.000 

7.02 

0.0375 

0.0065 

62.3 

120.0 

1.25250 

0.00 

0.000 
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run  # 

Wing 
radius  (m) 

Radius  of 
non-wash 
region  (m) 

Flapping 
Freq  (Hz) 

(j)max  - 
(t)min  (deg) 

Chord 

Factor 

yfactor 

Phase 
Shift:  |3 
(deg) 

Max. 

camber 

7.03 

0.0425 

0.0115 

62.3 

120.0 

1.25250 

0.00 

0.000 

7.04 

0.0375 

0.0065 

62.3 

150.0 

1.50000 

0.00 

0.000 

7.05 

0.0375 

0.0065 

62.3 

151.7 

1.00000 

1.25000 

-15.30 

0.000 

7.05.1 

0.0375 

0.0065 

62.3 

151.7 

1.25000 

1.25000 

-15.30 

0.000 

7.05.2 

0.0375 

0.0065 

62.3 

151.7 

1.50000 

1.25000 

-15.30 

0.000 

7.06 

0.0375 

0.0065 

62.3 

151.7 

1.00000 

1.25000 

-15.30 

0.000 

7.06.1 

0.0375 

0.0065 

62.3 

151.7 

1.25000 

1.25000 

-15.30 

0.000 

7.06.2 

0.0375 

0.0065 

62.3 

151.7 

1.50000 

1.25000 

-15.30 

0.000 

7.07 

0.0375 

0.0065 

62.3 

160.0 

1.50000 

1.60000 

0.00 

0.000 

7.08 

0.0375 

0.0065 

62.3 

151.7 

1.25000 

1.25000 

-15.30 

4.000 

7.09 

0.0375 

0.0065 

62.3 

151.7 

1.25000 

1.25000 

-15.30 

4.000 

prop 

0.0375 

0.0000 

145.8 

180.0 

1.21050 

7.12 

0.0375 

0.0000 

145.8 

180.0 

1.21050 

7.12.1 

0.0375 

0.0000 

145.8 

180.0 

1.21050 

7.10 

0.0375 

0.0000 

145.8 

180.0 

1.21050 

7.10 

0.0375 

0.0000 

145.8 

180.0 

1.21050 

7.10 

0.0375 

0.0000 

80.0 

180.0 

1.21050 

7.10 

0.0375 

0.0000 

80.0 

180.0 

1.21050 

7.13 

0.0375 

0.0000 

80.0 

160.0 

1.21050 

7.13 

0.0375 

0.0000 

80.0 

160.0 

1.21050 

7.14 

0.0375 

0.0000 

80.0 

160.0 

0.89193 

7.14 

0.0375 

0.0000 

80.0 

160.0 

0.89193 

7.15 

0.0375 

0.0000 

80.0 

180.0 

1.21050 

7.15 

0.0375 

0.0000 

80.0 

180.0 

1.21050 

7.16 

0.0375 

0.0000 

80.0 

140.0 

1.21050 

7.16 

0.0375 

0.0000 

80.0 

140.0 

1.21050 

7.17 

0.0375 

0.0000 

80.0 

180.0 

1.21050 

7.17 

0.0375 

0.0000 

80.0 

180.0 

1.21050 

7.18 

0.0375 

0.0000 

80.0 

180.0 

0.89193 

7.18 

0.0375 

0.0000 

80.0 

180.0 

0.89193 

7.19 

0.0375 

0.0000 

80.0 

160.0 

0.89193 

7.19 

0.0375 

0.0000 

80.0 

160.0 

0.89193 

7.20 

0.0375 

0.0000 

80.0 

160.0 

0.89193 

7.20 

0.0375 

0.0000 

80.0 

160.0 

0.89193 

7.20.1 

0.0375 

0.0000 

80.0 

160.0 

0.89193 

7.20.1 

0.0375 

0.0000 

80.0 

160.0 

0.89193 

7.21 

0.0375 

0.0000 

80.0 

160.0 

0.89193 

7.21 

0.0375 

0.0000 

80.0 

160.0 

0.89193 

dros/runl4 

0.0023 

0.0000 

211.3 

156.6 

left  wing 

0.000 

dros/runl4 

0.0023 

0.0000 

211.3 

157.8 

right  wing 

0.000 

dro  s/run  14 

0.0023 

0.0000 

211.3 

314.4 

both  wings 

0.000 
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Table  3.  Results  from  CFD  simulations 


Freq. 

required  for 
thrust  (w  two 

Mean  Thrust  wings)  equal  two 


run  # 

T  per  wing 
(N) 

Paero  Pef  wing 

(W) 

to  weight  - 

/,o(Hz) 

wings,  scaled 
to/i„ 

A 

''^swept 

(m^) 

Induced 
power  (W) 

Eigure  of 
merit 

1.01 

6.2600E-02 

l.OlOOE+00 

70.78 

1.3989 

0.00351 

0.3344 

23.91% 

1.1 

4.8100E-02 

5.9910E-01 

80.74 

1.2320 

0.00351 

0.3344 

27.15% 

1.2.2 

1.3098E-02 

1.1471E-01 

116.05 

1.6601 

0.00237 

0.4069 

24.51% 

1.2.3 

1.3533E-02 

1.1424E-01 

114.17 

1.5741 

0.00237 

0.4069 

25.85% 

2.01 

2.4900E-02 

2.6620E-01 

87.39 

1.4697 

0.00361 

0.3295 

22.42% 

2.02 

2.8400E-02 

3.4790E-01 

81.83 

1.5769 

0.00361 

0.3295 

20.90% 

2.03 

2.9400E-02 

4.5500E-01 

80.43 

1.9580 

0.00361 

0.3295 

16.83% 

2.04 

2.5300E-02 

4.1300E-01 

86.70 

2.2263 

0.00361 

0.3295 

14.80% 

2.05 

2.8000E-02 

4.9600E-01 

82.42 

2.2965 

0.00361 

0.3295 

14.35% 

2.06 

2.8700E-02 

81.40 

0.00361 

0.3295 

2.07 

2.0700E-02 

5.6070E-01 

95.85 

4.0841 

0.00361 

0.3295 

8.07% 

2.08 

2.2400E-02 

6.2380E-01 

92.14 

4.0364 

0.00361 

0.3295 

8.16% 

2.09 

2.1900E-02 

6.7710E-01 

93.19 

4.5322 

0.00361 

0.3295 

7.27% 

5.01 

3.9200E-02 

6.6140E-01 

69.65 

1.8487 

0.00361 

0.3295 

17.83% 

5.02 

2.9280E-02 

2.9700E-01 

80.59 

1.2859 

0.00361 

0.3295 

25.63% 

5.03 

3.1750E-02 

8.6430E-01 

77.40 

3.3142 

0.00361 

0.3295 

9.94% 

5.04 

3.6650E-02 

6.1610E-01 

72.04 

1.9049 

0.00361 

0.3295 

17.30% 

5.05 

3.1880E-02 

5.6010E-01 

77.24 

2.1346 

0.00361 

0.3295 

15.44% 

5.06 

3.2290E-02 

6.7840E-01 

76.75 

2.5363 

0.00361 

0.3295 

12.99% 

5.07 

3.4040E-02 

8.9080E-01 

74.75 

3.0769 

0.00361 

0.3295 

10.71% 

5.08 

3.8530E-02 

7.6630E-01 

70.26 

2.1980 

0.00361 

0.3295 

14.99% 

5.09 

3.8190E-02 

4.7440E-01 

70.57 

1.3789 

0.00361 

0.3295 

23.90% 

5.10 

3.9670E-02 

6.6130E-01 

69.24 

1.8156 

0.00361 

0.3295 

18.15% 

5.11 

3.2290E-02 

2.9940E-01 

76.75 

1.1194 

0.00361 

0.3295 

29.44% 

5.12 

3.5080E-02 

7.9960E-01 

73.63 

2.6400 

0.00361 

0.3295 

12.48% 

5.13 

3.7420E-02 

6.0460E-01 

71.29 

1.8119 

0.00361 

0.3295 

18.19% 

5.14 

3.7490E-02 

6.0070E-01 

71.22 

1.7952 

0.00361 

0.3295 

18.36% 

5.15 

3.9180E-02 

6.2830E-01 

69.67 

1.7575 

0.00361 

0.3295 

18.75% 

5.09.1 

3.3890E-02 

4.5570E-01 

74.91 

1.5845 

0.00361 

0.3295 

20.80% 

5.10.1 

3.5720E-02 

6.4720E-01 

72.97 

2.0797 

0.00361 

0.3295 

15.85% 

5.04.1 

4.9440E-02 

9.7450E-01 

72.18 

1.9230 

0.00361 

0.3295 

17.14% 

5.04.2 

3.2690E-02 

4.1770E-01 

76.27 

1.5331 

0.00361 

0.3295 

21.49% 

5.04.3 

4.8150E-02 

7.6150E-01 

62.85 

1.5635 

0.00361 

0.3295 

21.08% 

5.04.4 

7.9030E-02 

1.5548E+00 

49.06 

1.5181 

0.00525 

0.2733 

18.00% 

6.01 

2.2940E-02 

4.1130E-01 

91.05 

2.5680 

0.00361 

0.3295 

12.83% 

6.02 

1.8440E-02 

2.2700E-01 

101.56 

1.9666 

0.00361 

0.3295 

16.76% 

6.03 

1.1570E-02 

9.8000E-02 

128.21 

1.7082 

0.00361 

0.3295 

19.29% 

6.04 

6.5000E-03 

5.1700E-02 

171.05 

2.1402 

0.00361 

0.3295 

15.40% 

6.05 

2.1960E-02 

1.6830E-01 

93.06 

1.1219 

0.00361 

0.3295 

29.37% 

6.06 

2.4100E-02 

1.6490E-01 

88.83 

0.9561 

0.00361 

0.3295 

34.46% 

7.01 

1.4000E-02 

1.3460E-01 

116.55 

1.7627 

0.00286 

0.3705 

21.02% 

7.02 

1.5690E-02 

1.2750E-01 

110.10 

1.4073 

0.00286 

0.3705 

26.33% 

7.03 

2.2740E-02 

2.0950E-01 

91.45 

1.3253 

0.00351 

0.3344 

25.24% 

7.04 

3.1320E-02 

3.0890E-01 

77.92 

1.2090 

0.00357 

0.3314 

27.41% 

7.05 

2.1470E-02 

1.6350E-01 

94.12 

1.1274 

0.00361 

0.3295 

29.23% 
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run  # 

Mean  Thrust 

T  per  wing 
(N) 

Paero  Pef  wing 

(W) 

Ereq. 

required  for 
thrust  (w  two 
wings)  equal 
to  weight  - 

/,o(Hz) 

for  two 
wings,  scaled 
to/i„ 

A 

''^swept 

(m^) 

Induced 
power  (W) 

Eigure  of 
merit 

7.05.1 

2.6210E-02 

2.0900E-01 

85.18 

1.0685 

0.00361 

0.3295 

30.84% 

7.05.2 

3.1390E-02 

2.7270E-01 

77.84 

1.0637 

0.00361 

0.3295 

30.98% 

7.06 

2.2330E-02 

1.6160E-01 

92.29 

1.0506 

0.00361 

0.3295 

31.37% 

7.06.1 

2.6760E-02 

2.1920E-01 

84.30 

1.0863 

0.00361 

0.3295 

30.34% 

7.06.2 

3.2070E-02 

2.8990E-01 

77.01 

1.0950 

0.00361 

0.3295 

30.09% 

7.07 

1.4020E-02 

1.2940E-01 

116.47 

1.6910 

0.00381 

0.3209 

18.98% 

7.08 

2.8910E-02 

2.1900E-01 

81.11 

0.9665 

0.00361 

0.3295 

34.10% 

7.09 

3.0480E-02 

2.3000E-01 

78.99 

0.9376 

0.00361 

0.3295 

35.15% 

prop 

5.4085E-02 

3.2025E-01 

138.81 

0.5523 

0.00442 

0.2979 

53.94% 

7.12 

8.3805E-02 

5.0265E-01 

111.51 

0.4495 

0.00442 

0.2979 

66.29% 

7.12.1 

7.4799E-02 

4.7108E-01 

118.03 

0.4995 

0.00442 

0.2979 

59.64% 

7.10 

1.7099E-01 

6.2776E+00 

78.07 

1.9260 

0.00442 

0.2979 

15.47% 

7.10 

1.0698E-01 

1.0500E+00 

98.70 

0.6510 

0.00442 

0.2979 

45.77% 

7.10 

4.7248E-02 

9.7107E-01 

81.47 

2.0511 

0.00442 

0.2979 

14.53% 

7.10 

3.1514E-02 

1.5897E-01 

99.76 

0.6165 

0.00442 

0.2979 

48.33% 

7.13 

3.7690E-02 

7.8454E-01 

91.22 

2.3259 

0.00393 

0.3160 

13.59% 

7.13 

2.8644E-02 

1.3253E-01 

104.63 

0.5930 

0.00393 

0.3160 

53.29% 

7.14 

4.4053E-02 

5.6256E-01 

84.37 

1.3199 

0.00393 

0.3160 

23.94% 

7.14 

3.5150E-02 

1.4387E-01 

94.45 

0.4736 

0.00393 

0.3160 

66.73% 

7.15 

5.9089E-02 

1.7133E+00 

72.85 

2.5876 

0.00442 

0.2979 

11.51% 

7.15 

4.1714E-02 

2.2066E-01 

86.71 

0.5619 

0.00442 

0.2979 

53.03% 

7.16 

3.8161E-02 

1.8130E+00 

90.65 

5.2760 

0.00344 

0.3378 

6.40% 

7.16 

7.1886E-02 

1.3149E+00 

66.05 

1.4800 

0.00344 

0.3378 

22.83% 

7.17 

2.7796E-02 

3.6952E-01 

106.22 

1.7297 

0.00442 

0.2979 

17.22% 

7.17 

2.6579E-02 

1.3602E-01 

108.62 

0.6810 

0.00442 

0.2979 

43.75% 

7.18 

1.7579E-02 

4.0412E-01 

133.56 

3.7613 

0.00442 

0.2979 

7.92% 

7.18 

2.0368E-02 

1.1597E-01 

124.08 

0.8654 

0.00442 

0.2979 

34.43% 

7.19 

4.6718E-02 

7.1766E-01 

81.93 

1.5417 

0.00393 

0.3160 

20.50% 

7.19 

4.3185E-02 

3.8414E-01 

85.22 

0.9286 

0.00393 

0.3160 

34.03% 

7.20 

3.4265E-02 

3.4373E-01 

95.67 

1.1756 

0.00393 

0.3160 

26.88% 

7.20 

3.0294E-02 

1.3464E-01 

101.74 

0.5539 

0.00393 

0.3160 

57.05% 

7.20.1 

3.0555E-02 

3.1388E-01 

101.31 

1.2749 

0.00393 

0.3160 

24.79% 

7.20.1 

2.6597E-02 

1.1719E-01 

108.58 

0.5861 

0.00393 

0.3160 

53.92% 

7.21 

4.3312E-02 

4.7575E-01 

85.09 

1.1450 

0.00393 

0.3160 

27.60% 

7.21 

4.9317E-02 

3.5886E-01 

79.74 

0.7108 

0.00393 

0.3160 

44.46% 

1.44585E- 

dros/runl4 

6.4210E-06 

1.3355E-05 

05 

2.7621E-06 

20.68% 

1.45647E- 

dro  s/run  14 

6.4920E-06 

1.3543E-05 

05 

2.7978E-06 

20.66% 

2.90233E- 

dros/runl4 

1.2913E-05 

2.6898E-05 

05 

5.5598E-06 

20.67% 
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